/*
Copyright 2025 PerfXLab (Beijing) Technologies Co., Ltd.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
#include "OpenSparseBLAS_func.h"
#include <iostream>
#include <cstddef>
#include <complex>

// -----------------------------------------------------------------------------
// 通用稀疏 axpy 模板：y[indx[i]*incy] += alpha * x[i]
// -----------------------------------------------------------------------------
template <typename T>
void OpenSpB_axpy(
    int nz,
    const T& alpha,
    const T* x,
    const int* indx,
    T* y,
    int incy,
    OpenSpB_base_type index_base)
{
    // BLAS 风格的快速返回
    if (nz == 0 || alpha == T{}) return;

    // 合法性检查（仅当需要真正计算时）
    if (nz < 0 || incy <= 0) {
        std::cerr << "error: illegal input\n";
        return;
    }
    if (!(index_base == OPENSPBLAS_ZERO_BASE || index_base == OPENSPBLAS_ONE_BASE)) {
        std::cerr << "error: unsupported index_base\n";
        return;
    }
    if (!x || !indx || !y) {
        std::cerr << "error: null pointer\n";
        return;
    }

    const int base = (index_base == OPENSPBLAS_ONE_BASE) ? 1 : 0;

    for (int i = 0; i < nz; ++i) {
        const int idx = indx[i] - base;  // 允许 0/1 基
        y[idx * incy] += alpha * x[i];
    }
}

// -----------------------------------------------------------------------------
// 单精度实数 saxpy
// -----------------------------------------------------------------------------
void OpenSpB_saxpy(
    int nz,
    float alpha,
    const float* x, const int* indx,
    float* y, int incy,
    OpenSpB_base_type index_base)
{
    OpenSpB_axpy<float>(nz, alpha, x, indx, y, incy, index_base);
}

// -----------------------------------------------------------------------------
// 双精度实数 daxpy
// -----------------------------------------------------------------------------
void OpenSpB_daxpy(
    int nz,
    double alpha,
    const double* x, const int* indx,
    double* y, int incy,
    OpenSpB_base_type index_base)
{
    OpenSpB_axpy<double>(nz, alpha, x, indx, y, incy, index_base);
}

// -----------------------------------------------------------------------------
// 单精度复数 caxpy
// -----------------------------------------------------------------------------
void OpenSpB_caxpy(
    int nz,
    const void* alpha,
    const void* x, const int* indx,
    void* y, int incy,
    OpenSpB_base_type index_base)
{
    const auto a = *static_cast<const std::complex<float>*>(alpha);
    OpenSpB_axpy<std::complex<float>>(
        nz, a,
        static_cast<const std::complex<float>*>(x),
        indx,
        static_cast<std::complex<float>*>(y),
        incy,
        index_base);
}

// -----------------------------------------------------------------------------
// 双精度复数 zaxpy
// -----------------------------------------------------------------------------
void OpenSpB_zaxpy(
    int nz,
    const void* alpha,
    const void* x, const int* indx,
    void* y, int incy,
    OpenSpB_base_type index_base)
{
    const auto a = *static_cast<const std::complex<double>*>(alpha);
    OpenSpB_axpy<std::complex<double>>(
        nz, a,
        static_cast<const std::complex<double>*>(x),
        indx,
        static_cast<std::complex<double>*>(y),
        incy,
        index_base);
}
